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Abstract. Stochastic modeling has become a popular approach to quantify uncertainty in flows 
through heterogeneous porous media. The uncertainty in heterogeneous structure properties is of- 
ten parameterized by a high-dimensional random variable. This leads to a deterministic problem 
in a high-dimensional parameter space and the numerical computation becomes very challengeable 
as the dimension of the parameter space increases. To efficiently tackle the high-dimensionality, 
we propose a hybrid high dimensional model representation (HDMR) technique, through which the 
high-dimensional stochastic model is decomposed into a moderate-dimensional stochastic model in a 
most active random space and a few one-dimensional stochastic models. The derived low-dimensional 
stochastic models are solved by incorporating sparse grid stochastic collocation method into the pro- 
posed hybrid HDMR. The porous media properties such as permeability are often heterogeneous. 
To treat the heterogeneity, we use a mixed multiscale finite element method (MMsFEM) to simulate 
each of derived stochastic models. To capture the non-local spatial features of the porous media 
and the important effects of random variables, we can hierarchically incorporate the global informa- 
tion individually from each of random parameters. This significantly enhances the accuracy of the 
multiscale simulation. The synergy of the hybrid HDMR and the MMsFEM reduces the stochastic 
model of flows in both stochastic space and physical space, and significantly decreases the compu- 
tation complexity. We carefully analyze the proposed HDMR technique and the derived stochastic 
MMsFEM. A few numerical experiments are carried out for two-phase flows in random porous media 
and support the efficiency and accuracy of the MMsFEM based on the hybrid HDMR. 
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1. Introduction. The modeling of dynamic flow and transport processes in ge- 
ologic porous media plays a significant role in the management of natural resources 
such as in oil reservoirs and water aquifers. These porous media are often created 
by complex geological processes and may contain materials with widely varying abil- 
ities to transmit fluids. Thus multiscale phenomena occur in the model. Because of 
measurement errors and limited knowledge of physical properties, modeling of sub- 
surface flow and transport often uses random fields (e.g., permeability). The model's 
output can be accurately predicted by efficiently simulating the stochastic multiscale 
model. The existence of heterogeneity at multiple scales combined with this uncer- 
tainty brings significant challenges to the development of efficient algorithms for the 
simulation of the dynamic processes in random porous media. Therefore, the inter- 
est in developing stochastic multiscale methods for stochastic subsurface models has 
steadily grown in recent years [TOl |Tfl [14 l \19 \ \20 \ 126 ) [35] . 

The existence of uncertainty in random porous media is an important challenge 
for simulations. One way to describe the uncertainty is to model the random porous 
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media as a random field which satisfies certain statistical correlations. This naturally 
results in describing the flow and transport problem using stochastic partial differ- 
ential equations (SPDE). Many numerical methods have been developed for solving 
SPDE in the past decades. The existing stochastic numerical approaches roughly fall 
into the two categories: (1) non-intrusive schemes and (2) intrusive schemes. In non- 
intrusive schemes, the existing deterministic solvers are used without any modification 
to solve a (large) set of deterministic problems, which correspond to a set of samples 
from the random space. This leads to a set of outputs, which are used to recover 
the desired statistical quantities. Monte Carlo methods [5] and stochastic collocation 
methods fall into this category. In intrusive schemes, approximation of the probabilis- 
tic dependence is embedded directly in an existing deterministic method. Thus, the 
resulting numerical model is fundamentally different from the original deterministic 
model. Typical examples from this class arc stochastic Galcrkin |27j and perturbation 
methods [35]. A broad survey of these methods can be found in [35]. Among these 
methods, stochastic collocation methods have been paid much attention in the re- 
search community because they share the fast convergence property of the stochastic 
finite element methods while having the decoupled nature of Monte Carlo methods. 
A stochastic collocation method consists of two components: collocation points and 
interpolation operator. The collocation points are selected based on special rules. 
The different choice of collocation points leads to the different collocation methods 
including full-tensor product method [3] and Smolyak sparse grid method [4[ [29l [32] . 
Sparse grid stochastic collocation method is known to have the same asymptotic ac- 
curacy as full-tensor product collocation method by using fewer collocation points. 
However, the number of collocation points required in the sparse grid method will 
increase drastically for high-dimensional stochastic problems. Hence the method is 
still limited to the problems in a moderate-dimensional stochastic space. 

Due to the small correlation length and large variability of the covariance struc- 
ture, the uncertainties of random porous media property are usually parameterized 
by a large number of random variables, e.g., using truncated Fourier type expansion 
for random fields. Consequently, the model's input is defined in a high-dimensional 
random parameter space. Sampling in high-dimensional random space is very time- 
consuming. If the sampling of random space is conducted in full random space through 
stochastic collocation method, then the number of samples increases sharply with 
respect to the dimension of the random space. This is the notorious curse of di- 
mensionality, which poses great difficulties for the stochastic approximation in the 
high-dimensional stochastic space. Dimension reduction techniques such as high di- 
mensional model representation (HDMR) |31j , proper orthogonal decomposition and 
principal component analysis, are very helpful to alleviate the difficulty. Among these 
methods, HDMR is one of the efficient stochastic dimension reduction techniques and 
has caught extensive attention of researchers in recent years [25j [26l EH 033]- The 
HDMR was originally developed in the application of chemical models by [31] and 
now becomes a general set of quantitative assessment and analysis tools for capturing 
the high-dimensional relationship between model inputs and model outputs. HDMR 
has been used for improving the efficiency of deducing high-dimensional input-output 
system behavior, and can be employed to relieve the computation effort. Instead 
of dealing with the full high-dimensional random space, truncated HDMR technique 
can decompose a high-dimensional model into a set of low-dimensional models, each 
of which needs much less computational effort. The curse of dimensionality can be 
considerably suppressed by using HDMR. However, there exist some drawbacks in 
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traditional truncated HDMR [31] and adaptive HDMR gSJ For example, 

too many low-dimensional models need to be computed, and high-order of coop- 
erative effects from important random variables are neglected. To overcome these 
drawbacks, we propose a hybrid HDMR, which implicitly uses complete HDMR in 
a moderate-dimensional space and explicitly uses the firs-order truncated HDMR for 
the remaining dimensions. Thus the hybrid HDMR decomposes a high-dimensional 
model into a moderate-dimensional model and a few one-dimensional models. The 
moderate-dimensional space is spanned by the most active random dimensions. We 
use sensitivity analysis to identify the most active random dimensions. Then each 
low-dimensional stochastic model in the hybrid HDMR is solved by using sparse grid 
stochastic collocation method. The proposed hybrid HDMR renders a good approxi- 
mation in stochastic space and substantially improves efficiency for high-dimensional 
stochastic models. Therefore, a good trade-off between computation complexity and 
dimension reduction error is achieved in the hybrid HDMR technique. 

Porous media often have inherent complex heterogeneities and this poses a great 
challenge for numerical simulations. The heterogeneities may vary at very different 
scales. Simulating the flows in the porous media that use a very fine grid to cap- 
ture the heterogeneities is computationally very expensive and possibly infeasible. 
However, disregarding the heterogeneities can lead to large errors. To take into the 
heterogeneities, many multiscale methods (see [5J HJ [T7] and references there) have 
been developed in recent twenty years and successfully used to simulate the flows in 
heterogeneous porous media. Among these multiscale methods, MMsFEM [6j [T] has 
become an important numerical multiscale method. The main idea of MMsFEM is 
to incorporate the small-scale information into finite element velocity basis functions 
and capture their effect on the large scale via mixed finite element formulation. The 
MMsFEM retains local conservation of velocity flux and has been found to be use- 
ful for solving flow and transport equations in heterogenous porous media. In many 
cases, the multiscale basis functions can be computed overhead and used repeatedly 
in subsequent computations with different source terms, boundary conditions. This 
leads to a large computational saving in simulating the flow and transport process 
where the flow equation needs to be solved many times dynamically. When porous 
media exhibit non-separable scales, some global information is needed for representing 
non-local effects (e.g., channel, fracture and shale barriers) and used to construct mul- 
tiscale basis functions. Using global information can significantly improve accuracy 
and render agreement with geological realism p] [2j [18j [19] . If the global information 
is incorporated into MMsFEM, we refer to the MMsFEM as global MMsFEM (G- 
MMsFEM). In the paper, we extend the central concept of the global information to 
MMsFEM based on the hybrid HDMR. 

We incorporate hybrid HDMR, sparse grid stochastic collocation method and 
MMsFEM together and develop a stochastic dimension reduction MMsFEM. For the 
MMsFEM, we hierarchically utilize approximate global information, which is in tunc 
with the different parts of the hybrid HDMR. This implies the cheap computation for 
the approximate global information. The stochastic dimension reduction MMsFEM 
is able to reduce a high-dimensional stochastic multiscale model in both stochastic 
space and spatial physical space. Compared with traditional truncated HDMR tech- 
niques, much better efficiency and very good accuracy are achieved in the hybrid 
HDMR. We carefully analyze the proposed stochastic dimension reduction MMsFEM 
and investigate its applications for flows in random and heterogeneous porous media. 
Some important statistical properties (e.g., mean and variance) of the outputs of the 
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stochastic flow models are discussed. 

The rest of the paper is organized as follows. In Section 2, we briefly introduce the 
background on the flow and transport system in random porous media. In Section 3, 
we present a general framework of HDMR and propose the hybrid HDMR technique. 
Some theoretical results and computation complexity are also addressed in the section. 
Section 4 is devoted to presenting MMsFEM based on the hybrid HDMR. In Section 5, 
a few numerical examples are presented to demonstrate the performance of MMsFEM 
based on the hybrid HDMR. Finally, some conclusions and closing remarks are made. 

2. Background and notations. 

2.1. Two-phase flow system and its stochastic parametrization. Let D 

be a convex bounded domain in R d (d = 2,3) and (f2, F, P) be a probability space, 
where il is the set of outcomes, F is the u-algebra generated by f2, and P is a 
probability measure. 

Let k(x,uj) be a random permeability field. We consider two-phase flow and 
transportation in the random permeability field. Here the two phases is referred as 
water and oil, designated by subscripts w and o, respectively. The equations of two- 
phase flow and transport can be given (in the absence of gravity and capillary effects) 
by flow equations 

- div(X(S)k(x, io)Vp) = q, (2.1) 

— + div(uf w (S)) = 0, (2.2) 
at 

where the total mobility X(S) is given by X(S) = X W (S) + X (S) and q is a source term. 
Here X W (S) — k rw {S)/ \i w and X (S) = k ro (S)/ /j , where \i and \i w are viscosities of 
oil and water phases, respectively, and k rw (S) and k ro (S) are relative permeabilities 
of oil and water phases, respectively. Here f w (S) is the fractional flow of water and 
given by /„, = X W /(X W + A D ), The equation (|2.ip is the flow equation and the equation 
(12.21) is the transport (or saturation) equation. According to Darcy's law, the total 
velocity u in (|2.2p is given by 

u = u w + u a = —X(S)kX7p. (2-3) 

A random field can be usually parameterized by Fourier type expansion such as 
expansion of Karhunen-Loeve (ref.[23]), polynomial chaos and wavelet [22] . This often 
gives rise to an infinite-dimensional random space. For computation, we truncate such 
an expansion to approximate the random field. Then k(x, w) can be formally described 

by 

k{x,uj) « k{x,6i(u),6 2 (w),...,6 N (u)). (2.4) 

For example, if a random field is characterized by a covariance structure, then the 
random field can be approximated by a finite sum of uncorrelatcd random variables 
through truncated Karhunen-Loeve expansion (KLE). To get accurate approxima- 
tion, a large number N of random parameters in (|2.4[) arc required. This leads to a 
deterministic model in a high-dimensional random parameter space. For simplicity of 
presentation, we make the following assumption in the paper, 

k(x, uj) = k(x, 0), where 8 := 0(w) = (0(wi), 6 2 {ui), 6 N {u)) e R N . 

Let I N C R N be the image of 6, i.e., I N = and p(6) = n^iM#i) 

be the joint probability function of (6±, ■ ■ ■ ,6n). By equations (|2.1[) . (|2.2[) . (|2.3p and 
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parametrization of permeability field, we formulate the following stochastic two-phase 
flow system: find random fields p(x, 9, t), v(x, 6, t), S(x, Q,t) : D x I N x (0, T] — > K 
such that they almost surely (a.s) satisfy the following equations subject to initial 
condition and boundary condition 

div(ii) = q 

u = -X(S)k(x,Q)V P (2.5) 
— + div(uf w (S))=0. 

For simplicity, we will suppress the spatial variable x and temporal variable t in 
the rest of paper when no ambiguity occurs. 

2.2. Sparse grid stochastic collocation method. For stochastic two-phase 
flow systems (|2.5p . the statistic properties (e.g., mean and variance) of solutions are 
of the quantities of interest. These properties may be obtained by first sampling 
the parameter random space using, for example, a Monte Carlo method or sparse 
grid collocation method, then solving the deterministic problems for the samples and 
analyzing the corresponding results to obtain the desired statistic quantities. The 
convergence of Monte Carlo methods is slow and a very large number of samples may 
be required, which leads to high computational cost. Instead, we use Smolyak sparse 
grid collocation method [35], where the Smolyak interpolant A(N + £, N) (£ > 1) is a 
linear combination of tensor product interpolations with the property: only products 
with a relatively small number of nodes are used and the linear combination is chosen 
in such a way that an interpolation property for N = 1 is preserved for TV > 1 [3]. 
In the notation A(N + £,N), the £ represents the interpolation level. Sparse grid 
collocation method is known to have the same asymptotic accuracy as tensor product 
collocation method, while requiring many fewer collocation points as the parameter 
dimension increases [29j . 

Sparse grids have been successfully applied to stochastic collocation in many re- 
cent works (e.g., [25, 26, 19,). Based on Smolyak formula (ref.[4]), a set of collocation 
points {G^y^fj in I N are specially chosen, where 7V C is the number of collocation 

points. With these chosen collocation points and corresponding weig hts {w^}f^, 
the statistical properties of the solutions can be obtained. At each of the colloca- 
tion points, the deterministic system (|2.5[) is solved and the output, for example, 
S(x,QV\t) is obtained. Then the mean of S(x, Q,t) can be estimated by 

E[s(x,e,t)]= I s{x,e,t) P (e)de 

J IN 

« / A(N + £,N)[S(x,e,t)]p(G)dG = V S(x,e U) ,t)w u) . 
JlN 



Here the weights {w^> are determined by the basis functions of A(N + £, N) and 
joint probability function p(O) (ref. [13]). Similarly, the variance of S(x,Q^,t) can 
be obtained by 

Var[S(x,e,t)]= [ (S(x,e,t)-E[S(x,e,t)]) 2 p(Q)dQ 



n c (2.6) 

s 2 (x, e u \t)w U) - e 2 [s(x, e, t)}. 

3 = 1 
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Let H(N + £, N) denote the number of collocation points for Smolyak sparse grid 
interpolation A(N + £,N). Then it follows that (see @]) 

2 e 

H{N + £,N) w — A £ foriV»l. (2.7) 

This implies that the number of collocation points algebraically increases with respect 
to the dimension N. We utilize Smolyak sparse grid collocation for the numerical 
computation. The stochastic approximation of the Smolyak sparse grid collocation 
method depends on the number of collocation points and the dimension N of the 
random parameter space. The convergence analysis in |29| implies that the conver- 
gence of Smolyak sparse grid collocation is exponential with respect to the number 
of Smolyak points, but depends on the parameter dimension N. This exponential 
convergence rate behaves algebraically for N 3> 1. 

3. High dimensional model representation. The truncated KLE leads to 
the two-phase flow deterministic model (|2.5|l with a high-dimensional parameter 0. 
The most challenge part of solving such a high-dimensional stochastic system is to 
discretize the high-dimensional random parameter. There exist a few methods for the 
discretization of the random space [HEHD^i 33 35 . Among these methods, sparse grid 
collocation method has been widely used and generates completely decoupled systems, 
each of which is the same size as the deterministic system. The sparse grid collocation 
method is usually very efficient in moderate-dimensional spaces. However, when the 
dimension of the random parameter is large and the approximation of stochastic space 
is conducted in the high-dimensional stochastic space, a large number of collocation 
points are required by (|2.7j) and we need to solve deterministic model (|2.5j) at all of 
the collocation points. Consequently, the high dimensionality will deteriorate the effi- 
ciency for the collocation method. To overcome the difficulty, we use high dimensional 
model representation (HDMR) to reduce the stochastic dimension and enhance the 
efficiency of simulation. By truncating HDMR, the high-dimensional model can be 
decomposed into a set of low-dimensional models. Hence the computation efforts can 
be significantly relieved [3SJ |33] . 

In the section, we present a general HDMR framework and propose a hybrid 
HDMR. 

3.1. A general HDMR framework. In the section, we adopt the concept on 
decompositions of multivariate functions in |21j to present a general HDMR framework 
in terms of operator theory. Let T be a linear space of real-valued functions /(0) 
defined on a cube I N and := (0\, 62, ■ • • , 0n) 6 I N ■ In the paper, /(0) represents 
a relationship between the input random vector and an output (e.g., saturation, 
water-cut). 

To introduce HDMR, we define a set of projection operators as follows. 
Definition 3.1. J21f Let {Pj}jLi be a set of commuting projection operators on 
T satisfying the following property: 

Pj(/) = / if f does not depend on 0j, and Pj(/) does not depend on 9j. (3.1) 
Let u C {1, • ■ ■ , N} be a subset and / the identity operator, we define 



P u :=n jeu Pj, and P := /. 
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The property (|3.1[) implies that Pi ... jv(/) is a constant. The identity operator / can 
be decomposed through the above commuting projections Pj (J — 1, • ■ • , N): 

I = Ilf =1 [P j + (I-P j )} = J2 ( U * n l iy )(ll,. v: a P 3 ). (3.2) 
Then the equation (|3.2p leads to a decomposition of / given by 

/= Yu /«> whcrc /" := ( n jeu( / - p j))(n ie{ i,...,jv}\uPj)/. (3-3) 

uC{l,-,JV} ^ ' ^ ' 

where / u depends only on variables with indices in u and f$ is a constant. The 
equation (|3.3p gives an abstract HDMR formulation. 

Associated with projection P u , we define P u := P{i,...,jv}\ u . Then P u / only 
depends on the variables with indices in u. As a result, the component / u of / 
defined in (|3.3j) can be recursively computed and explicitly computed, respectively by 
(ref. [U) 

/»=Pu/-£/v and / u = ^(-l) |u| - |v, Pv/. (3.4) 

vCu vCu 

Next we specify the projection operator in (|3.ip . Let \i be a measure on Borel 
subsets of I N and a product measure with unit mass, i.e., 



= 1, i = !,•••, N. (3.5) 



We assume that the functions in T are integrable with respect to \x. For any / € T 
and j e {1, • • • , TV}, we define 

P i f(Q) = Jf(Q)diM j (0 i ). (3.6) 

We can check that the projections in (|3.6p satisfy the conditions in (|3.ip . Conse- 
quently, for any u C {1, • • • , N} 



p u /(e)= / /(e)n^ u d M 

In particular, 



P /= / f(Q)dm- 

The inner product (•, ■) on J- induced by the measure fi is defined as follows: 



(f,h) := / f(e)h(Q)d»(0), f l9 eF. 

The norm || • \\ T on J 7 is defined by \\f\\ T := (/, f) 1 ' 2 . 

Definition 3.2. Let F$ and F u ($ ^ u C {1, ■ • • , N} ) be defined as follows: 

J-0:= {/ G T : f is constant C}, 

J- u := {/ S J- : f = f u depends on the variables with indices in u, and f £ nj eu ifer(Pj)}. 
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By the definition of the projection operator, we can obtain the following theorem, 
which gives a decomposition of J- . 

Theorem 3.3. fffffl Let J" u be defined in Definition \3.S\ for any u C {1, • ■ • , N}. 
Then 

(1) J 7 — ©uC{l,--- ,JV}-^u- 

(2) T n := (ll jeu (I -P^P^. 

For any u C {1, • • • , N}, we define 

M u = ^eu^-PjO^Pu. 

Then we can show the operators M u are commutative projection operators and mu- 
tually orthogonal, i.e., 

M u = M u , M U M V = for u ^ v. 
Theorem 13.31 gives the HDMR expansion of / by 

/= ]T f u = Yl M "/' where M u / e J- u . (3.7) 

uC{l,---,AT} uC{l,— ,N} 

Theorem 13.31 also implies that 

Wffr = E ll M »/l&> P0(M U /) = O foru^0. 

uC{l,- ,N} 

Any set of commutative projectors {M v } generate a distributive lattice whose ele- 
ments are obtained by all possible combinations (addition and multiplication) of the 
projectors in the set. The lattice has a unique maximal projection operator Ai, which 
gives the algebraically best approximation to the functions in J- |15j . The range of 
the maximal project operator Ai for the lattice is the union of the ranges of {M v }. 
Because the commutative projectors {M v } are mutually orthogonal here, the maxi- 
mal project operator M. and the range !Fm of M have their explicit expressions as 
follows: 

M = E M v, Tm = fflv^v 

V 

As more orthogonal projectors are retained in the set, the resultant approximation 
by its maximal projection operator M. will become better. 

By using the notion of maximal projection operator, we define the order of HDMR 
approximation. Consider the set of projectors {M u : u C {1, • • • , N}, |u| < q}, the 
best approximation of / is given by 

/ « M q f := J2 M u /- 

uc{i,— A ; 

u| <q 

This gives the q-th order HDMR approximation of /. General speaking, higher or- 
der HDMR approximations are never worse than lower order HDMR approximation. 
Adding new orthogonal projectors into a sum of orthogonal projectors always produce 
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a new maximal projection operator with better accuracy. It is argued in |31j that a 
low order HDMR (e.g., up to the third-order) usually gives a good approximation for 
most practical problems. 

If the measure fj. in (|3.5|) is taken to be the probability measure 

dfi(e) = P {&)de = u^piie^dOi, 

then the resulting HDMR defined in ([317]) is called ANOVA-HDMR [3T] . Fixed a 
point B := (01,02, " j0iv)- If the measure /i in (|3.5j) is taken as the Dirac measure 
located at the point 0, i.e., 

duQB) = n£ 1 5(0 i - 6i)d6 h 

then the resulting HDMR is Cut-HDMR. The point B is so called cut-point or anchor 
point. In ANOVA-HDMR, the P v f involves N — |v| dimensional integration and 
the computation of the components of / is expensive. In Cut-HDMR, the P v / only 
involve the evaluation in a N — |v| dimensional space and the computation is cheap 
and straightforward. Because of this reason, we use Cut-HDMR for the analysis and 
computation in the paper. 

Remark 3.1. Using \3J\ and induction method, it can be shown J21f that 

M q f = Y^ a u E ? »/> * rea « ■ = {-^) q ' j ( N ~^~ 3 \ (3-8) 

3=0 uC{l,-,JV} \ 1 3 J 

u | < q 

This gives a straightforward calculation for A4 q f . 

3.2. A hybrid HDMR. In this subsection, we follow the idea of abstract 
HDMR and present HDMR in a straightforward way. For practical computation, 
a hybrid HDMR is proposed for approximation of HDMR. 

By (|3.3p . the HDMR of /(0) can be represented in the form 

n 

/(B) = fo + J2 + fa( 9 *> W + ■■■ /ia-n(»i, fla, • • • , On)- (3-9) 

i— 1 i<j 

Here /o is the zeroth-order component denoting the mean effect of /(B). The first- 
order component /i(0j) represents the individual contribution of the input 6i and the 
second-order component fij(9i, 6j) represents the cooperative effects of 9i and 0j and 
so on. For most realistic physical systems, low-order HDMR M q f (e.g., q < 3) may 
give a good approximation [3T] . The g-th order truncated HDMR can be rewritten as 

q 

m — 1 i\<---i m 

The AWg/ can approximate /(6) through truncating the HDMR. The Ai q f consists 
of a large number of component terms for high-dimensional models, the computation 
for all the components may be costly. For many cases, some high-order cooperative 
effects can not be neglected in computation models, especially when a model strongly 
relies on a few dependent variables. Consequently, the traditional truncated HDMR 
in (|3.10l) may not give very good approximation. To eliminate or alleviate the above 
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drawbacks of traditional truncated HDMR, wc propose a new truncated HDMR, which 
is termed by hybrid HDMR. 

We use Fourier amplitude sensitivity test [5] to pick up the most active dimensions 
from the component of the high-dimensional random vector 0. Let {fj}jLi be the 
first-order components defined in Eq. (|3.9[) and a 2 denote variance. We may assume 
that <r 2 (fi) > <J 2 {f2) > ••• > c 2 (/aO- Otherwise, we re-order the index set {j}jLi 
by their monotonically decreasing values of variance. Alternatively, we can use their 
monotonically decreasing values of expectation, E[fi] > E[f2\ > •• ■ > E\fx\. In the 
paper, we restrict ourself to the variance sensitivity. We set a threshold constant with 
< C < 1 and find an optimal J such that 

£a 2 (^-)/f> 2 (^-)>C. (3.H) 

3=1 3=1 

Then we define to be the most J active dimensions. We note that a 2 (fj) gives 

us information on the impact of 9j when it is acting alone on the output. It is clear 
that if the change of the random input 9j within its range lead to a significant change 
on the output, then <J 2 (fj) is large. Therefore, the criterion (|3.11[) is reasonable to 
identify the most active dimensions. According to practical situation, we can adjust 
the value of £ in Eq. (|3.11[) such that J is much less than N. Wc define an index set 
J = {1, 2, • • ■ , J}. Then we define the hybrid HDMR as following, 

Mjf:=^M v f+ ( 3 - 12 ) 

vcj ie{i,- ,n}\j 

Here the set of projectors {M v ,v C J} U |Mj,i E {1, ■ ■ ■ 7 N} \ generate a 

lattice and its maximal project operator is M. j. By (|3.12[) . the hybrid HDMR Aijf 
consists of two parts: the first part X) v c j M v / is the complete HDMR on the most 
active dimensions indexed in Jj, the second part N}\j ^»/ ^ s ^ ne first-order 

truncated HDMR on the remaining dimensions {1, ••• ,N} \ J. By p.4p . we can 
rewrite equation (|3.12p by 

Mjf:=Pjf+ ( 3 - 13 ) 

»e{i,- ,n}\j 

We note that the operator P j projects the N- variable function / to a function defined 
on the J the most active dimensions. The term P jf gives the dominant contribution 
to M jf. Any first-order components are usually important, these components are re- 
tained in M. jf. In Cut-HDMR, there is no error for the approximation Ai jf of /(0) 
whenever the point is located the J-th dimensional subvolume {0 j, ■ • ■ , On} 
across the cut-point 0. 

If we use traditional truncated HDMR to approximate the term Pjf in (|3.13|) . 
then we can obtain the adaptive HDMR developed in (25j [34] , 

Mfj:= M v(Pj/)+ E /<(*«)• (3-14) 

v<ZJ,\v\<q ie{l,— ,N}\J 

For simplicity of notation, we will suppress q in M.°j q in the paper when the truncation 
order q is not emphasized. The following proposition gives the relationship among /, 
M jf and M a /f. 
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PROPOSITION 3.4. Let M jf and M a jf be defined in \3.1S\) and {3.1$ , respec- 
tively. Then 

||/ - M a j d f\\% =\\.f- Mjf\\l + \\Mjf Mff\\%. (3.15) 

Proof. We note the equality 

/ - Mff = (/ - Mjf) + (Mjf - Mff). (3.16) 

Because the set of HDMR components in / — M jf do not intersect with the set of 
HDMR components in Mjf - M a j f , Theorem |3~31 implies that 

(f-Mjf,M J f-M J d f)=0. (3.17) 

Then the equation (|3.15p follows immediately by combining equation (|3.16p and (|3 . 1 7[) . 
□ 

The equation (|3. 15|) implies that hybrid HDMR has better approximation than 
adaptive HDMR. 

The mean of / can be computed by summing the mean of all HDMR components 
of / for both ANOVA-HDMR and Cut-HDMR. The variance of / is the sum of the 
variance of HDMR components of / for ANOVA-HDMR. However, the direct sum- 
mation of variance of components of / may not equal to variance of / for Cut-HDMR. 
Consequently, we may want to have a truncated ANOVA-HDMR to approximate the 
variance of / for Cut-HDMR. If we directly derive a truncated ANOVA-HDMR, the 
computation is very costly and more expensive than the direct computation of the 
variance of / itself. This is because a truncated ANOVA-HDMR requires comput- 
ing many high-dimensional integrations. To overcome the difficulty, we can use an 
efficient two-step approach to derive a hybrid ANOVA-HDMR through the hybrid 
Cut-HDMR. 

Let O j be the J-dim variable with indices in J. Using the hybrid HDMR for- 
mulation (j3~T3l> . the hybrid Cut-HDMR f cut {9) of / has the following form 

/»*(0) = r\Qj) + E fi Ut &)' where/™* (Qj):=Pjf. (3.18) 

»e{i,- ,N}\J 

We use Mj to act on f cut (&) to get a hybrid ANOVA-HDMR / a "™ a (6>) of /, which 
has the following form 

f anova (Q):=Mj-f cut = f anova (Qj)+ fi nova (0i)- (3-19) 

i£{l,- ,N}\J 

Then we have the following theorem. 

Theorem 3.5. Let / CMt (6) and / a »°™(9) be defined in £3~T$\) and £3~W\) , re- 
spectively. Then 

E[f cut ] = E[f anova ] = E[f cut ] + E U! U % (3-20) 

ie{i,- ,n}\j 

Var[f cut ] = Var[.f anova ] = Var[f cut ] + M/fl- (3-21) 

ie{i,--- ,n}\j 
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Proof. For the proof, we need to calculate the HDMR components in J anoua (Q). 
By definition, we derive the term f anova (Qj') as follows, 

/- (&j) = ! Pjf cut = I r^ik^d^) 

JlN-J 

= f j cut (Q J )Ti lt jd l , l (e l )+ e [ /r'cw^^ 

/JY J »e{i,- ,n}\j JiN 

= r t {Q J )+ e 

-N}\J 

We define / anouQ = E[f cut \. Then for any i 6 {1, • ■ • , N} \ J, we have 
f~(0 t ) = [ f cut (0)n? =1 dfi s (6 s ) - f~ 

JlN-l s^i 

= f f cut (ej)iis=idvs(o s )+ V I /; ut (^)nf =1 ^ s (0 s )-/ o ~ 

je{i,-,JV}V7 (-323) 

= /r* ) + E[t ut ] + e E tf? ut ] - f~ 

je{i,- ,n}\j 

:= /f^i) + Ci. 

It is obvious that E[ff nova ] = by (pT2"3l) . Consequently, it follows that by (|3~2"2"j) . 
This proves the equality For AN OVA- HDMR, we note that 



Var[ < /* an01,a ] || yotiowo ||2 ||^anova||2 _j_ ^ ^ ||yanoua| 

ie{i,--- 

= Var[/ a " OUQ ] + V^[fr° Va ] 

iG{l,--- ,iV}\,7 

= Var[/ c " t ]+ Var [/r 4 ], 

»e{i,- 



(3.24) 



where equations (|3.22p and (|3.23p have been used in the last step. Further, by equa- 
tions (j3~!?2]) and (j3~2"3"]) we have 

/«*(e) = r\Qj) + e /<"'(*«) 

ie{i,--- -N}\J 

=r nova {Qj)- E ^[/r 4 ]+ E 

,N}\J ,N}\J 

= r nova (e)-( e (^[/r']+^; 

This implies that Var[/ cut ] = Var[/ ano ™]. Hence the proof is completed. □ 
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Theorem 13.51 implies that we can compute the mean and variance of f cut by 
directly summing the means and variances of the components of f cut . This can sig- 
nificantly reduce the complexity of computation. In addition, the derived hybrid 
ANOVA-HDMR f anova is easily obtain by using the two-step approach through the 
the hybrid Cut-HDMR f cut . Our further calculation shows that the traditional trun- 
cated Cut-HDMR (|3T0|) and adaptive Cut-HDMR (f3~T4| do not have these merits. 
We use sparse grid quadrature |13j to compute the mean and variance in the paper. 

Compared with the traditional truncated HDMR and adaptive HDMR, the hy- 
brid HDMR defined in Eq. ([3~T3"l) has the following advantages: (a) The hybrid HDMR 
consists of significantly fewer terms than the traditional truncated HDMR and adap- 
tive HDMR, the total computational effort of the hybrid HDMR can be substantially 
reduced, (b) Since the hybrid HDMR inherently includes the all cooperative contri- 
bution from the most active dimensions, approximation accuracy is not worse (maybe 
better) in hybrid HDMR than the traditional truncated HDMR and adaptive HDMR; 
(c) According to Theorem 13.51 the computation of variance of hybrid HDMR is much 
more efficient. 

3.3. Analysis of computation complexity. In the previous subsection, we 
have addressed the accuracy of the hybrid HDMR and made some comparisons with 
traditional truncated HDMR and adaptive HDMR. In this subsection, we discuss the 
computation efficiency for the various HDMR techniques when Smolyak sparse grid 
collocation is used. We remind that the HDMR means Cut-HDMR. 

Let C(/, £) be the number of sparse grid collocation points with level £ in full 
random parameter dimension space and C(M q f,£) the total number of sparse grid 
collocation points with level I in the traditional truncated HDMR Ai q f. We define 
C(M jf, t) and C(M a j q f, £) in the similar way. We first consider the case 1 = 2 and 
q = 2 and calculate the total number of sparse grid collocation points for the different 
approaches. By using (|2.7j) and (|3.8|) . we have for N ^> 1 

C(/,2) 
C(M 2 /,2) 

C(Mjf,2) 
C(M a j d 2 f,2) 

Consequently, if J ; 

C{Mjf,2) < C(f,2) < C(M a j d 2 f,2) < C(M 2 f,2). 

This means that the hybrid HDMR Ai jf requires the smallest number of collocation 
points when £ = 2 and q = 2. This case is our particular interest in the paper. 

Next we consider some other cases. By some subtle calculation, it can be shown 
that if J w 30 and N > (f then 

C(Mjf, 3) < C(M a j d J,3) < C(M2f,3) < C(/,3). 

This implies that the computation effort in the hybrid HDMR Mjf is the least 
for £ = 3 when the number of the most active dimensions is moderate for a high- 
dimensional problem. However, if the number J of the most active dimensions is 



= H(N + 2,N) rj 27V 2 , 

= £ (*) h u + 2 >j) = 5N + 13 ( 2 ) * 13iv2 ' 

= H(J + 2,J) + (N— J)H(1 + 2, 1) w 2 J 2 + 5(N - J), 

= £ ( J ) H U + 2 > 3) + ( N - J ) F ( 1 + 2 ' !) ~ 13 j2 + 5N - 

j= i v ? / 

s N/2, it follows immediately that 



14 



L. Jiang, J. D. Moulton and J. Wei 



large (e.g., J 3> 30), we can show that 

C(M a j d q f, t) < C(M q f, I) < C(Mjf, t) < C(f, £), where £ > 3 and q < 3. (3.26) 

The relationship p.26|) tells us that adaptive HDMR J\A°j f may be the most efficient 
if higher level collocation is required and J is large as well. 

For a high-dimensional stochastic model, if the number J of the most active 
dimensions is large and high level (e.g., level 3 and above) sparse grid collocation is 
required to approximate the term Pjf in (|3 . 1 3|) . using the adaptive HDMR (|3.14l) 
may improve efficiency by the comparison p.26[) . The adaptive HDMR has been 
extensively considered in many recent papers |16[ 1251 [2"6l 134) . By our experience, we 
find that the number of the most active dimensions J is often less than ^ as N is 

large. If Pjf is smooth with respect to Qj, low level sparse grid collocation (e.g., 
level 2) can usually provide accurate approximation. Many problems can fall into this 
situation and they are our special interest in the paper. 

3.4. Integrating HDMR and sparse grid collocation method. As stated 
in Subsection l2.2[ the sparse grid stochastic collocation method can reduce the stochas- 
tic two-phase flow system (|2.5p into a set of deterministic two-phase flow systems. 
However, it suffers from curse of dimensionality with increasing the number of the 
stochastic dimensions. Integrating HDMR and sparse grid collocation method is a 
practical way to overcome this difficulty and can significantly enhance the efficiency. 

Without loss of generality, we use saturation solution S(&) as an example to 
present the technique for the hybrid Cut-HDMR. By using (|3.13j) . (|3.4p and Smolyak 
sparse grid interpolation, we have 

S(0) « MjS(Q) = S(Oj) + Yl s ^ 

ie{i,...,N}\j 

= s(Oj)+ J2 m) - (n - j)s (3 27) 



ie{i,...,N}\j 
A(J + t,J)[S(Qj)] 4 



Y, Aii + t,i)[8(ei)]-(N-j)s , 

ie{i,...,N}\J 



where S(Qj) = PjS(Q), S(0i) = P l S(&) and S = P S , (6) = S(@). Here 6 is the 
cut-point for a Cut-HDMR. The choice of cut-point may effect on the approximation 
of truncated Cut-HDMR. The study in [12] argued that an optimal choice of the cut- 
point is the center point of a sparse grid quadrature. In the paper, we will use such a 
cut-point for computation. Due to Theorem I3.5[ the mean and variance of S can be 
approximated by 



E[S(Q)] « E 



A(J + £, J)[S(Qj)] 



(N-J)S , 



Var[5(0)] « Var 



A(J + £,J)[S(Qj)] 



ie{i,...,N}\j 



A(l+£,l)[S(6i)] 



(3.28) 



. : i n}\j 



A(l+l,l)[S(6i)] 



From (|3.27[) . we find that the approximation error of the mean and variance comes 
from the two sources: truncated HDMR and Smolyak sparse grid quadrature. Let 
{Qj,w^y^ii be the set of pairs of collocation points and weights in the most active 



MMsFEM based on hybrid HDMR 



15 



dimension space I J and {9\ \\ w^}^L 1 the set of pairs of collocation points and 
weights in /. We note that Nj = H{J + £, J) and N. L = H(l + £, 1). Then by (|cT2"gj) . 
the mean of S can be computed by 

Nj Ni 

E[S{<3)} » S(Of)w^ + J2 J2 S(e^ k) )w (k) — (N — J)S . (3.29) 

3=1 ie{l,...,N}\J k=l 

Similarly, we can compute Var[5(0)] in terms of evaluations of S(Qj) and S(8i) at 
the collocation points. 

Because each of terms in adaptive HDMR M°j q S are usually correlated each 
other, the variance of j\A a j q S is not equal to the summation of the variance of all 

terms in (jOij) . Let {Q^ ,wW}*=i (where N c = H(N+£,N)) be the set of pairs of 
collocation points and weights in full dimensions I . To compute the variance using 
adaptive HDMR, we need project each collocation point 0^ onto the components 
{Q [ v n) : v C J, |v| < q} and {0^ : i e {1, • • ■ , N} \ J} and interpolate all terms in 
(I3.14j) . Then we use (|2.6|) to calculate the variance. This process involves at most 
H (N + e,N)[N-J + £« =1 ( 3 7 )] Smolyak sparse grid interpolations. The computa- 
tion of these interpolations is usually very expensive when N and J are large. The 
computation of variance using adaptive HDMR is usually much more expensive than 
using hybrid HDMR. The numerical experiments in Section [5] will confirm the point. 

4. Mixed multiscale finite element method. In Section 13.41 we have dis- 
cussed integrating hybrid HDMR and sparse grid collocation method to reduce the 
computation complexity from high-dimensional stochastic spaces. The permeability 
field is often heterogeneous in porous media. It is necessary to use a numerical method 
to capture the heterogeneity. MMsFEM is one of such numerical methods and has 
been widely used in simulating flows in heterogeneous porous media [2j 1 1 9 j . To sim- 
ulate the two-phase flow system (|2.5[) , it is necessary to retain local conservation for 
velocity (or flux). To this end, we use MMsFEM to solve the flow equation and ob- 
tain locally conservative velocity. Using MMsFEM coarsens the multiscale model in 
spatial space and can significantly enhance the computation efficiency. 

Corresponding to the hybrid HDMR expansion of S(<d) in (|3.27[) . i.e., 

MjS(Q) = S(Qj) + J2 S(6i) - (N - J)S , 

ie{l,...,N}\J 

the velocity u(8) in (12.51) also admits the same hybrid HDMR expansion as following 
Mju(&) = u(Qj)+ Yl u(9i) - (N - J)u , (4.1) 

ie{l,...,N}\J 

where u(Qj) = Pju(@), u(8i) — Piu(6) and ito = u(Q). We use j) to obtain 
5(0j), u{9i) to obtain S{6i) and uq to obtain So- 

Without loss of generality, we may assume that the boundary condition in the 
flow equation of (|3.27[) is no flow boundary condition. Let k(Qj) — P jk(Q) and 
k(9i) = Pj/s(0). Then we can uniformly formulate the mixed formulation of equations 
of u{Qj) and u(Qi) as following, 

f (X(S)ky 1 u + Vp = inL> 

< div(u) = q in D (4.2) 
u ■ n = on dD. 
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Here ii(Qj) and u(6i) are corresponding to the coefficients k(Qj) and k(6i), respec- 
tively. Let fco = fc(0). Then uq is the solution to equation (|4.2p if k and S are replaced 
by ko and So, respectively. 

The weak mixed formulation of (|4.2[) reads: seek (u,p) G Ho(div,D) x L 2 (D)/R 
such that they satisfy the equation 

- (dw(v),p) = Vu G H (diw, £>) 
(g,r) VreL 2 {D). 

Let C i?o(div,D) and C L 2 (D)/R be the finite element spaces for velocity 
and pressure, respectively. Then the numerical mixed formulation of (|4.2j) is to find 
(uh,Ph) G Vh x Q/j such that they satisfy 



(4.3) 




(X(S)k) 1 Uh,Vhj - (div(v h ),Ph) = Vvh G Vh 
(dxv(u h ),r h ) = (q,r h ) Mr h 6 Qh- 



We use MMsFEM for (|4.3[) . It means that mixed finite element approximation 
is performed on coarse grid, where the multiscale basis functions arc defined. In 
MMsFEM, we use piccewise constant basis functions on coarse grid for pressure. For 
the velocity, we define multiscale velocity basis functions. The degree of freedom of 
multiscale velocity basis function is defined on interface of coarse grid. Let ef be a 
generic edge or face of the coarse block K. The velocity multiscale basis equation 
associated with ef is defined by 

-div(kVwf ) = mK 



-kVw, 



K . n _ } J a K bf nds 



- on ef (4.4) 



else. 



For local mixed MsFEM [6], bf = n, the normal vector. If the media demonstrate 
strong non-local features including channels, fracture and shale barriers, some limited 
global information is needed to define the boundary condition bf to improve accuracy 
of approximation [2 [18] . We will specify the boundary condition bf for different parts 
in (|4.ip . Then ipf = —kS/wf defines multiscale velocity basis function associated to 
ef , and the multiscale finite dimensional space for velocity is defined by 



K,i 



K 



It is well-known that using approximate single-phase global velocity information 
can considerably improve accuracy for multiscale simulation of two-phase flows [TJ [2] 
[TS1QI5]. For the two-phase flow system (|2.5p . the single-phase global velocity u sg (Q) 
solves the following equation 

(fc(9))~V S9 + V Ps9 = inD 

div(it sg ) = q in D (4.5) 
Usg ■ n = on dD. 
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By using hybrid HDMR and sparse grid interpolation, the u sg (Q) admits the following 
approximation 

Mju sg {Q) = u sg {Qj)+ Y u sg {0i) - {N - J)u sgfi 

ie{i,...,N}\j 

~ <u sg>0 + Y Usg,j(0j) > + Y ii sg (9 t ) - (N - J)u sgfi 
^ j£J ' ie{l,-,N}\J 

= \0--J)u S g,o + Y {is a( e o)\ + Y u sg {6i)-{N- J)u sgfi 
^ jej > ie{i,...,N}\j 

(i - j)u sg , + Y AO- + *, i) ) 

+ Y A(l+£',l)[u sg (6i)]-{N-J)u sg , 
ie{i,...,N}\j 

where ii sg (9j) = Pjii sg (<d) and u sg (9j) = PiU sg (Q). Here the interpolation levels £ 
and £' can be different. Because {9j}j£j are the most active dimensions, it is usually 
desirable that I > £' . Now we are ready to specifically describe the multiscale finite 
element space for u(Q j), u(9i) (i £ {1, . . . , N} \ J) and u . 

• To construct the multiscale basis functions for u(Q j), we take k = k(Q j) in 

and 

bf = b?{Qj)= ({l-J)u ag » + YAl+^)[u sg {e 3 )])\ e K. (4.6) 

The multiscale finite element space for u(Q j) is defined by 

Vh(®j) = (B<P?(ej), where ^(Qj) = -%{®j)Vw?{®j). (4.7) 

K,i 

• To construct the multiscale basis functions for u{9i) (i € {1, . . . , N} we 
take k = k(9i) in P~4]l and 

bf = bf(9,) = ^(1+Al)[fl s5 (^)])| ef . (4.8) 
The multiscale finite element space for u(9i) is defined by 

Vhtfi) = iff (ft), where (ft) = -fc(ft)Vu;f (ft). (4.9) 



K,i 

• To construct the multiscale basis functions for uq, we replace k by fc(O) in 
(|4~4]) and 

h = u sg,o\ef ■ 

The multiscale finite element space for uq is defined by 

Vh(&) = V>f (6), where (6) = -fc(S)Vt«f (6). 
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For different parts of the hybrid HDMR expansion M.ju(&) defined in (|4.1|) , 
we use different boundary conditions for multiscale basis equations. This increases 
the hierarchies of multiscale basis functions, which arc in tunc with the sensitivity of 
random parameter dimensions. 

By (|3.29l) , to compute the moments of outputs we only need to construct the mul- 
tiscale finite element space Vh(Qj) at the set of collocation point {6j }^ J j, where 
Nj = H(J + £, J). For arbitrary Qj £ I J , the boundary condition bf(Qj) is com- 
pletely determined by u sg {9^), where j G J and m = 1, . . . , H(l + £, 1). Here {9^} 
is a set of collocation points in I. Similarly, we also need to construct the multiscale 
finite element space Vh(9i) at the set of 1-dimensional collocation points {o'f ■*}, where 
i € {1, . . . , N} \ J and j' = 1, . . . , H(l + £, 1). For arbitrary 9i G /, the boundary con- 
dition bf(9i) is completely determined by u sg (9^ n ^), where m' = 1, . . . , H(l + £', 1) 
and {9j ^} is a set of collocation points in /. These single-phase global velocity in- 
formation {u sg (9j m ^), u sg {9 t i n ■*)} is pre-computed and can be repeatedly used in the 
whole stochastic flow simulations. 

The convergence analysis for the mixed MsFEM using approximate global infor- 
mation can be found in jTSJ Q]5] . From the above discussion, here the approximate 
global information is completely determined by stochastic single-phase velocity con- 
tributed from individual random dimensions 9i (i = 1,...,N). According to the 
analysis in [18], the resonance error between the coarse mesh size h and the effect of 
individual contributions of 9i {i = 1, . . . , N) is negligible. This error is a major source 
using local MMsFEM |19| . However, the resonance error between the coarse mesh 
size h and the effect of cooperative contributions of is retained, but is usually small. 
This is our motivation to use limited global information determined in 1-dimensional 
random space. If the random permeability field has low variance and does not exhibit 
global features, we can use local mixed MsFEM for the stochastic simulation. 

Define 

u h (x, 0) := A( J+£, J) [u h (x, Qj)) + J2 A( - 1+ ^ ^ e *)] -(N-J)u , h ^)- 

i£{l,...,N}\J 

For velocity, we actually compute Uh(x, 0), namely, the numerical solutions u^(x, O j) 
and Uh(x,6i) arc evaluated at their collocation points. Wc define u{x, 0) in the same 
way. Now we analyze the total error between u(x, 0) and Uh(x, 0). By triangle 
inequality, 

E[\\u(x,Q)-u h (x,Q)\\ LHD) ] 

< E[\\u(x, 0) - Mju(x, 0)|U= P) ] + E[\\Mju(x, 0) - u(x, Q)\\ LHD) ] 
+ E[\\u(x,G) -u h (x,e)\\ L 2 (D) ] 

• £hdmr ~\~ &col ~l~ £msi 

where Shdmr represents the error from hybrid HDMR, £ co i the error by sparse grid 
collocation and £ ms the error of MMsFEM. The error Etdmr depends on the choice 
of the most active dimensions J ', the error £ co i depends on the number of collocation 
points and random dimensions, and the error £ ms relies on the coarse mesh size h 
and the approximation of the boundary conditions (|4.6[) and (|4.8p to u sg (@j) and 
u sg (9i), respectively. According to the works in [HUH], if the boundary conditions in 
(14. 61) and (|4.8|) are replaced by u sg (Qj) and u sg (9i), respectively, then the error £ ms 
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only depends on coarse mesh size h. The error £ co i has been extensively discussed in 
many recent literatures, e.g., [3J|31|33], and we focus on the errors Shdmr and £ ms m 
the paper. 

5. Numerical results. In this section, we present a few numerical results for 
two-phase flow in random porous media. The hybrid HDMR and MMsFEM are 
incorporated together to enhance simulation efficiency. If the random porous media 
have some non-local spatial features and high variance, we use limited global informa- 
tion to construct multiscale velocity information to achieve better accuracy. We will 
combine MMsFEM and different truncated HDMR techniques (hybrid HDMR M j 
and adaptive HDMR Aij 1 ) for the flow simulations. The accuracy and efficiency of 
computation will be compared for the different methods. 

For numerical simulation, we assume that k(x 1 0) to be a logarithmic random 
field, i.e., k(x, 0) := exp(a(x, 0)). This assumption assures that k(x, O) is positive 
and the well-posedness of flow equation in (|2.5|) . Here a(x, uj) is a stochastic field and 
can be characterized by its covariance function cov[a]: D x D — > R by 

cov[o](a;i,X2) ■= E[(a(xi) - E[a(x 1 )])(a(x 2 ) - E[a(x 2 )})] . 

For the numerical experiments, we use a two point exponential covariance function 
for the stochastic field a, 

r u \ 2 f \xi~x 2 \ 2 | y\ - y 2 | 2 \ , K ^ 

cov[a\(x 1 ,y 1 ;x 2 ,y 2 ) = a exp ^ — — J, (5.1) 

where (xi, yi) (i = 1, 2) is the spatial coordinate in 2D, a 1 is the variance of stochastic 
field a, and l x and l y denote the correlation length in the x— and y— direction, re- 
spectively. We will specify these parameters for the covariance function in numerical 
examples. Using Karhuncn-Loeve expansion (KLE) the random field a(x, 0) admits 
the following decomposition 

N 

a(x, 0) := E[a] + J2 V^M^iM, (5-2) 

i=l 

where the random vector := (0i,6 2 ,--- ,0n) G and the random variables 
{di{ijj)}^ =1 are mutually orthogonal and have zero mean and unit variance. The 
effects of permeability field a(x, 0) with uniform, beta and Gaussian distributions on 
mean and standard deviation of output were discussed in |23) . where the numerical 
results showed the similar peak values of standard deviation for the three different 
distributions. Therefore, in the section we assume that each 9i is i.i.d. uniform on 
[— 1 , 1] . This will not lose the main feature of the output uncertainty and hurt the 
performance of the proposed methods. 

When MMsFEM is used, the fine grid is coarsened to form a uniform coarse grid. 
We solve the pressure equation on the coarse grid using MMsFEM and then recon- 
struct the fine-scale velocity field as a superposition of the multiscale basis functions. 
The reconstructed velocity field is used to solve the saturation equation with a fi- 
nite volume method on the fine grid. We solve the two-phase flow system (|2.5p by 
the classical IMPES (Implicit Pressure Explicit Saturation). The temporal discreti- 
sation is an implicit scheme, which is unconditionally stable but leads to a nonlinear 
system (Newton- Raphson iteration solves the nonlinear system). In all numerical sim- 
ulations, mixed multiscale basis functions are constructed once at the beginning of 
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computation. In the discussion, we refer to the grid where multiscale basis functions 
are constructed as the coarse grid. The limited global information is computed on the 
fine grid. 

We compare the saturation fields and water-cut data as a function of pore volume 
injected (PVT). The water-cut is defined as the fraction of water in the produced fluid 
and is given by q w /qt, where qt = q +qw, with q and q w being the flow rates of oil and 
water at the production edge of the model. In particular, q w = J dDou t f(S)u-nds, qt = 
JgQout u ■ nds, where dD out is the out-flow boundary. Pore volume injected is defined 
as PVI = y- J Q qt(r)dT, with V p being the total pore volume of the system, provides 
the dimensionless time for the displacement. We consider a traditional quarter five- 
spot problem (e.g., pQ) on a square domain D, where the water is injected at left top 
corner and oil is produced at the right lower corner of the rectangular domain (in the 
2D examples). 

5.1. Linear transportation. In this subsection, we consider the case when the 
transport (saturation) equation is linear in the flow system (|2.5p . For this purpose, 
we consider the mobility of water and mobility of oil defined as following, 

X W (S) = S, \ (S) = l-S. 

Consequently, the total mobility X(S) = 1 and the fractional flow of water f w {S) = S. 
Then the two-phase flow system (|2.5p reduces to a linear flow model. Since X(S) = 1, 
the velocity is not updated to compute saturation. We use the example to investigate 
the performance of the hybrid HDMR technique and MMsFEM for linear models. 

The stochastic permeability field k(x, O) is given by k(x,Q) = exp (a(x, 6)) , 
where a(x, O) is characterized by the covariance function cov[a] defined in (|5.ip . whose 
parameters is defined as following: a 1 = 1, l x = l y = 0.1. We truncate the KLE (|5.2|) 
after the first 80 terms to represent the random field a(x, 0) with E[a] = 1. This 
implies that k(x, O) is defined in an 80-dimcnsional random parameter space, i.e., 
N = 80. Fig. 15.11 depicts a realization of the random field a(x, O). The stochastic 
field k(x, O) is defined in a 60 x 60 fine grid. We choose 6x6 coarse grid and apply 
MMsFEM to compute velocity. The time step is taken to be 0.02 PVI for discretizing 
temporal variable. 

To make comparison, we use various methods to simulate the linear two-phase 
flow system: fine-scale mixed finite element method (MFEM) on full random dimen- 
sions, MMsFEM on full random dimensions, MMsFEM based on adaptive HDMR, 
and MMsFEM based on hybrid HDMR. The solution computed by fine-scale MFEM 
is referred to as the reference solution in the paper. To assess the performance of 
MMsFEM and HDMR, we compute the mean and standard deviation (std) for the 
quantifies of interest from the two-phase flow model, such as saturation and water- 
cut. Since the random field in the example does not have strong non-local features, 
the local MMsFEM (L-MMsFEM) generally gives an accurate approximation [18]. 
Hence the local MMsFEM is utilized for simulation in this example. To choose the 
most active random dimensions for adaptive HDMR and hybrid HDMR, we choose a 
threshold constant £ = 0.9 and use the criterion (|3.11[) . This produces 31 most active 
dimensions. Smolyak sparse grid collocation method with level 2 is used to tackle 
stochastic space. We evaluate the flow model's outputs at the collocation points and 
use the associate weights to compute the mean and standard deviation of the outputs. 
Table 15.11 lists the number of deterministic models to be solved for the above four dif- 
ferent methods. From the table, we find that MFEM on full random dimensions needs 
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to compute the largest number of deterministic models (fine-scale models) , MMsFEM 
based on hybrid HDMR needs to compute the smallest number of deterministic mod- 
els (coarse-scale models). This means that MMsFEM based on hybrid HDMR can 
significantly reduce the computation complexity. 



Table 5.1 

number of deterministic models for the different methods 



methods 


number of deterministic models 


MFEM on full random dim. 


12961 fine-scale models 


MMsFEM on full random dim 


12961 coarse-scale models 


MMsFEM based on adaptive HDMR 


6446 coarse-scale models 


MMsFEM based on hybrid HDMR 


2231 coarse-scale models 



Fig. 15.21 shows the point-wise mean of saturation map at PVI = 0.4. From the 
figure, we observe that both L-MMsFEM based on hybrid HDMR and L-MMsFEM 
based on adaptive HDMR provide very accurate approximation for the mean of sat- 
uration. They are almost to identical to the reference mean computed by fine-scale 
MMFEM on full random dimensions. Fig. 15.31 illustrates the point-wise standard 
deviation of saturation at PVI = 0.4. By the figure, we have four observations: (1) 
Compared with reference standard deviation, the approximations of standard devia- 
tion by the three multiscale models (L-MMsFEM on full random dim., L-MMsFEM 
based on hybrid HDMR and L-MMsFEM based on adaptive HDMR) are good; (2) 
L-MMsFEM on full random dim. renders better approximation for standard de- 
viation than L-MMsFEM based on truncated HDMR techniques; (3) L-MMsFEM 
based on hybrid HDMR gives almost the same standard deviation as the L-MMsFEM 
based on adaptive HDMR; (4) The variance mostly occurs around the front of water 
flow. Fig. 15.41 depicts the mean of water-cut curves. We see that the four water-cut 
curves are almost identical. This demonstrates that a good approximation has been 
achieved. Fig. 15.51 depicts the standard deviation of water-cut curves by the four dif- 
ferent methods. By the figure, we see that both L-MMsFEM based on hybrid HDMR 
and L-MMsFEM based on adaptive HDMR can render a very good approximation for 
standard deviation at most time instances. The standard deviation of water-cut curve 
by L-MMsFEM based on hybrid HDMR is almost identical to the standard deviation 
of water-cut curve by L-MMsFEM based on adaptive HDMR. 

Next we discuss the relative errors of mean and standard deviation in the sense of 
integration with regarding to physical domain. Let S r /W r , S ms j/W ms j, S mSiaH /W ms , a H 
and Sms,hH /Wms,hH be the saturation/water-cut using MFEM on full random dimen- 
sions, MMsFEM on full random dimensions, MMsFEM based on adaptive HDMR and 
MMsFEM based hybrid HDMR, respectively. Then we define the relative errors of 
saturation mean and saturation standard deviation between S r and S ms .f as following, 

gm , s) = \\E[S r }(x) - E[S msJ }(x)\\ LHD) = \\a[S r ](x)-cr[S m s,f](x)\\ LHD) 

\\ E [Sr](x)\\ L i(D) h[Sr](x)\\ L i( D ) 

(5.3) 

where a denotes the standard deviation operation. We can similarly define the relative 
errors £™ H {S) and £^{S) between S r and S mSt aH, and the relative errors f^(S') 
and £fffl(S) between S r and S ms ,hH- We list these errors on saturation at 0.4 PVI in 
Table O From the table, we see that L-MMsFEM based on hybrid HDMR gives the 
same accuracy of mean as L-MMsFEM based on adaptive HDMR. They are both very 
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A realization of random field a(x, 0) 




Fig. 5.1. A realization of random field a(x, ©) 

close to the error of mean by L-MMsFEM on full random dimensions. The relative 
error of standard deviation by L-MMsFEM based on hybrid HDMR also has a good 
agreement with that of L-MMsFEM based on adaptive HDMR. The table implies that 
the errors only from adaptive HDMR and hybrid HDMR are very small. The main 
source of errors is from spatial multiscalc approximation. 



Table 5.2 

relative errors of mean and std. on saturation (at 0.4 PVI) 



methods 


relative error of mean 


relative error of std. 


L-MMsFEM on full random dim. 


1.223540e-002 


5.302182e-002 


L-MMsFEM based on adaptive HDMR 


1.361910c-002 


9.389291c-002 


L-MMsFEM based on hybrid HDMR 


1.361910e-002 


9.482789e-002 



For water-cut, the relative errors of mean and standard deviation between W r 
and W ms j are defined by 

fm (w] = \\m r ](t) - E[W ms , f ](t)\\L*ao,i] fstd(m = \\a[W r }(t)-a[W msJ }(t)\\ LHl0A]) 

ms[ 1 \\E[W r ](t)\\ L2{m) ' mA ' h[W r ](t)\\ LH[0A]) 

(5.4) 

where [0, 1] is the temporal domain. The relative errors of water-cut such as £"^ L H (W), 
£^h{W), £™ h {W) and £^(W), arc defined in a similar way. We list these relative 
errors of water-cut in Table I5T51 The table shows the similar behavior of these methods 
to the situation of saturation. The performance of L-MMsFEM based on hybrid 
HDMR is almost the same as the performance of L-MMsFEM based on adaptive 
HDMR. 



Table 5.3 

relative errors of mean and std. on water- cut 



methods 


relative error of mean 


relative error of std. 


L-MMsFEM on full random dim. 


1.445791e-002 


4.339328e-002 


L-MMsFEM based on adaptive HDMR 


1.856908c-002 


1.006136e-001 


L-MMsFEM based on hybrid HDMR 


1.856908c-002 


1.040452e-001 
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fine-scale MFEM on full random dim. L-MMsFEM on full random dim. 




Fig. 5.2. Point-wise mean of saturation map at PVI=0.4- Top left: MFEM on full random 
dimensions; Top right: L-MMsFEM on full random dimensions; Bottom left: L-MMsFEM based on 
hybrid HDMR; Bottom right: L-MMsFEM based on adaptive HDMR 



The number of the most active dimensions may have an important impact on 
computation efficiency and accuracy. To this end, we choose some different threshold 
constants in (pTTTj) . Specifically, we take C = 0.75, 0.8, 0.85, 0.9, 0.95. These threshold 
constants correspond to the number of the most active dimensions dim(Q j) = 1 9, 
22, 26, 31, 41. Fig. [O] illustrates the CUP times of L-MMsFEM based on adap- 
tive HDMR and L-MMsFEM based on hybrid HDMR for the different number of the 
most active dimensions. From the figure, we get two observations: (1) CPU time 
of L-MMsFEM based on hybrid HDMR is only a fraction of the CUP time using 
L-MMsFEM based on adaptive HDMR; (2) As the number of the most active di- 
mensions increases, the CPU time of L-MMsFEM based on hybrid HDMR increases 
mildly, but the CPU time of L-MMsFEM based on adaptive HDMR increases dras- 
tically. By the comparison, we see that hybrid HDMR is much more efficient than 
adaptive HDMR. There are two reasons for adaptive HDMR to take more CPU time: 
(1) adaptive HDMR needs to solve more deterministic models; (2) to compute stan- 
dard deviation (or variance), a large number of stochastic interpolations are involved. 
Fig. 15.71 depicts the relative errors of mean (left) and standard deviation (right) for 
saturation at 0.4 PV1. By the figure, we find that increasing the number of the most 
active dimensions can substantially improve the accuracy of the mean and standard 
deviation for saturation. We also see that the L-MMsFEM based on hybrid HDMR 
has almost the same accuracy as the L-MMsFEM based on adaptive HDMR for sat- 
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Fig. 5.3. Point-wise standard deviation of saturation contour at PVI=0.4- Top left: MFEM 
on full random dimensions; Top right: L-MMsFEM on full random dimensions; Bottom left: L- 
MMsFEM based on hybrid HDMR; Bottom right: L-MMsFEM based on adaptive HDMR 



uration approximation. Fig. 15.81 describes the the relative errors of water-cut mean 
(left) and water-cut standard deviation (right) for the different number of the most 
active dimensions. From the figure, we observe that the effect of increasing dim(Qj) 
may not be as sensitive as the case of saturation for enhancement of accuracy. The 
reason is that the MMsFEM error is the main contribution of total error and has 
greater impact on the water-cut than HDMR itself. 

5.2. Non-linear transportation. In the subsection, we consider the flow sys- 
tem (|2.5[) where the saturation is nonlinear. In the example, the mobility of water 
and the mobility of oil arc defined by 

X w (S) = S 2 / f i w , A„(5) = (I - Sf/no. 

Here \i w j 1 \i Q = 0.1, the ration between viscosity of water and oil. Consequently, the 
fractional function f w (S) of water is given by 

S 2 

fw{S) = S 2 + 0.1(1 -S) 2 " 

This results in a non-linear two-phase flow system. 

We again consider the random permeability k(x, 0) = exp (a(x, 0)). Here co- 
variance function cov[a] of a(x, 0) is defined in ()5.1[) with l x = l y = 0.2 and cr = 1. 
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PVI time 



Fig. 5.4. Mean of water-cut curves 




PVI time 

Fig. 5.5. Standard deviation of water- cut curves 

The mean of a(x, 0) is highly heterogenous and its map is depicted in Fig. 15.91 (left). 
It is actually obtained by extracting and rescaling an SPE 10 [7] permeability field 
(the 45-th layer) . We truncate the KLE (|5.2[) after the first 30 terms to represent the 
random field a(x, 0) and so the random field k(x, 0) is defined in an 30-dimensional 
random parameter space. Fig. 15.91 depicts a realization of the random field a(x, 0) 
(right). The stochastic field k(x, 0) is defined in a 60 x 60 fine grid. We choose 6x6 
coarse grid for MMsFEM to compute velocity. To discretize the temporal variable of 
saturation equation, the time step is taken to be 0.01 PVI. 

Fig. 15.91 shows that the permeability field k(x, 0) exhibit some channel features, 
which have important impacts on flows simulation. To achieve accurate simulation re- 
sults, limited global information can improve the accuracy [T|l2lll8|. In the subsection 
we incorporate hybrid HDMR technique and adaptive HDMR technique with both 
local MMsFEM (L-MMsFEM) and global MMsFEM (G-MMsFEM) and develop dif- 
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CUP time vs. number of most active dimensions 




Fig. 5.6. CPU times for the different number of the most active dimensions, 19, 22, 26, 31, 41. 
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Fig. 5.7. Relative errors of mean (left) and standard deviation (right) for saturation at 0.4 PVI. 

ferent methods. We will make some comparisons for these methods and explore their 
performance. To get the most active random dimensions for the truncated HDMR 
tcchnqiues, we choose a threshold constant £ = 0.9 and use the criterion (|3.11|) . This 
gives rise to 14 most active dimensions in the 30-dim random parameter space. We 
use the collocation points and weights associated with level 2 Smolyak sparse grid 
collocation method to compute mean and standard deviation. 

Fig. 15.101 shows the point- wise mean of saturations at 0.4 PVI for different meth- 
ods. We can clearly see that G-MMsFEM on full random dimensions has the best 
approximation to the reference solution given by fine-scale MFEM on full random 
dimensions. The figure also demonstrates that the multiscale methods with the trun- 
cated HDMR techniques provide good approximation. Fig. 15.111 describes the point- 
wise standard deviation of saturation at 0.4 PVI for those different methods. From 
the figure, we get two observations: (1) the large variance occurs around of front of 
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Fig. 5.8. Relative errors of mean (left) and standard deviation (right) for water-cut. 




Fig. 5.9. A realization of random field a(x, 0) 



water flow; (2) the hybrid HDMR technique gives almost the same standard devia- 
tion map as the adaptive HDMR technique. We note that these behaviors share the 
similarities in Subsection 15. II We also compute the water-cut curves for the different 
methods. For a better visualization, here we only present the water-cut curves by 
the four methods: MFEM on full random dim., G-MMsFEM on full random dim., 
G-MMsFEM based on adaptive HDMR and G-MMsFEM based on hybrid HDMR. 
Fig. I5.12l dcpicts the mean of water-cut curves for the four methods. From the figure, 
we see that the four curves overlap each other at almost all times. We note that there 
exists some fluctuation right after water break-through time. The reason may be that 
the value of water-cut changes sharply right after water break-through time. The 
standard deviation of the water-cut curves for the four methods are illustrated in Fig. 
15.131 By the figure, we find that the variance of water-cut mostly occurs between the 
water-break-through time and time when water is full of the reservoir. 

In order to carefully measure the differences caused by L-MMsFEM and G- 
MMsFEM integrating with adative HDMR and hybrid HDMR, we follow the con- 
cept in Subsection 15.11 and compute the relative errors between the reference solution 
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Fine-scale MFEM on full random dim. G-MMsFEM based on hybrid HDMR G-MMsFEM based on adaptive HDMR 
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Fig. 5.10. Point-wise mean of saturations at 0.4 PVI for different methods 



Fine-scale MFEM on full random dim. G-MMsFEM based on hybrid HDMR G-MMsFEM based on adaptive HDMR 
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Fig. 5.11. Point-wise standard deviation of saturations at 0.4 PVI for different methods 



and the solutions of the various coarse models (MMsFEM with (or without) HDMR 
technqiues) . Here the reference solution is solved by fine-scale MFEM on full random 
dimensions. The relative errors of saturation and water-cut are defined in a similar 
way as in (|5.3[) and (|5.3[) , respectively. Table 15.41 and 15.51 list the relative errors of 
mean and standard deviation for saturation at 0.4 PVI and water-cut, respectively. 
From the two tables, we find: (1) limited global information can enhance the accuracy 
of MMsFEM; (2) the results by adaptive and hybrid HDMR are very close to each 
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Fig. 5.12. Mean of water-cut curves for MFEM on full random dim., G-MMsFEM on full 
random dim., G-MMsFEM based on adaptive HDMR and G-MMsFEM based on hybrid HDMR 




Fig. 5.13. Standard deviation of water-cut curves for MFEM on full random dim., G-MMsFEM 
on full random dim., G-MMsFEM based on adaptive HDMR and G-MMsFEM based on hybrid 
HDMR 



other. 

Table 5.4 

relative errors of mean and std. on saturation (at 0-4 PVI) 



methods 


relative error of mean 


relative error of std. 


L-MMsFEM on full random dim. 


4.170726e-002 


1.545560c-001 


G-MMsFEM on full random dim. 


8.950902e-003 


4.027147e-002 


L-MMsFEM based on adaptive HDMR 


4.788991c-002 


2.051337e-001 


G-MMsFEM based on adaptive HDMR 


3.283281e-002 


1.608348e-001 


L-MMsFEM based on hybrid HDMR 


4.788991e-002 


2.013145e-001 


G-MMsFEM based on hybrid HDMR 


3.283281e-002 


1.520236c-001 



Finally we discuss the efficiency for the different approaches. To the end, we 
record the CPU times for the seven different methods: fine-scale MFEM on full 
random dimensions, L-MMsFEM on full random dimensions, G-MMsFEM on full 
random dimensions, L-MMsFEM based on adaptive HDMR, G-MMsFEM based on 
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Table 5.5 

relative errors of mean and std. on water- cut 



methods 


relative error of mean 


relative error of std. 


L-MMsFEM on full random dim. 


1.320466e-002 


8.836354e-002 


G-MMsFEM on full random dim. 


8.000854e-003 


4.298053e-002 


L-MMsFEM based on adaptive HDMR 


1.423242e-002 


1.352706e-001 


G-MMsFEM based on adaptive HDMR 


1.162734e-002 


1.212131e-001 


L-MMsFEM based on hybrid HDMR 


1.423242e-002 


1.447519e-001 


G-MMsFEM based on hybrid HDMR 


1.162734e-002 


1.248755c-001 




Fig. 5.14. CPU times for different approximation models 

adaptive HDMR, L-MMsFEM based on hybrid HDMR and G-MMsFEM based on 
hybrid HDMR. Fig. 15.141 shows the CUP times for the seven different approaches. 
We see that the approaches of hybrid HDMR need the least CUP time and achieves 
the best efficiency. The CPU time of hybrid HDMR is only a fraction of the CUP 
time of adaptive HDMR. In the meanwhile, the CUP time of G-MMsFEM is slightly 
larger than the CUP time of L-MMsFEM when truncated HDMR techniques are syn- 
thesized. For adaptive HDMR, a lot of CUP time is spent in computing variance. 
However, computation of variance in hybrid HDMR is straightforward and requires 
very little CPU time. This shows the merit of the hybrid HDMR. 

6. Conclusions. In the paper, we presented a framework of high-dimensional 
model representation (HDMR) for stochastic MMsFEM. A hybrid HDMR technique 
has been proposed and it decomposes a high-dimensional stochastic model into a 
moderate-dimensional stochastic model and a few one-dimensional stochastic models. 
Some criteria were presented to determine the moderate dimensions. We have devel- 
oped MMsFEM based on hybrid HDMR by integrating it with MMsFEM. MMsFEM 
based on hybrid HDMR can significantly reduce the original model's complexity in 
both multiscale physical space and high-dimensional stochastic space. In the frame- 
work, we also presented MMsFEM based on adaptive HDMR, which has been widely 
used in stochastic model reduction. Compared with adaptive HDMR, the hybrid 
HDMR is much more efficient and retain the the same (or better) accuracy as adap- 
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tive HDMR. The computation time of the hybrid HDMR is a fraction of the com- 
putation time of adaptive HDMR. To capture strong non-local features in multiscalc 
models, we have incorporated important global information into multiscale computa- 
tion. This can improve the approximation accuracy of the proposed coarse multiscale 
models. We carefully analyzed the proposed MMsFEM using HDMR techniques and 
discussed computation efficiency and approximation errors. The MMsFEM based on 
HDMR techniques was applied to flows in random porous media. Both linear and 
non-linear stochastic flows models were discussed and the numerical results confirmed 
the performance of the proposed approaches. 
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